Goldstein-price function
Contributed by: Benoît Legat
In this example, we consider the minimization of the Goldstein-price function.
using SumOfSquares
using DynamicPolynomialsCreate symbolic variables (not JuMP decision variables)
@polyvar x[1:2](DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}[x₁, x₂],)To use Sum-of-Squares Programming, we first need to pick an SDP solver, see here for a list of the available choices.
import Clarabel
using Dualization
model = SOSModel(dual_optimizer(Clarabel.Optimizer))A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Dual model with Clarabel attachedCreate a JuMP decision variable for the lower bound
@variable(model, γ)\[ γ \]
f(x) is the Goldstein-Price function
f1 = x[1] + x[2] + 1
f2 = 19 - 14*x[1] + 3*x[1]^2 - 14*x[2] + 6*x[1]*x[2] + 3*x[2]^2
f3 = 2*x[1] - 3*x[2]
f4 = 18 - 32*x[1] + 12*x[1]^2 + 48*x[2] - 36*x[1]*x[2] + 27*x[2]^2
f = (1 + f1^2*f2) * (30 + f3^2*f4)\[ 600 + 720x_{2} + 720x_{1} + 3060x_{2}^{2} - 4680x_{1}x_{2} + 1260x_{1}^{2} + 12288x_{2}^{3} - 19296x_{1}x_{2}^{2} + 7344x_{1}^{2}x_{2} - 1072x_{1}^{3} + 14346x_{2}^{4} - 23616x_{1}x_{2}^{3} + 7776x_{1}^{2}x_{2}^{2} + 5784x_{1}^{3}x_{2} - 2454x_{1}^{4} + 1944x_{2}^{5} - 11880x_{1}x_{2}^{4} + 5040x_{1}^{2}x_{2}^{3} + 9840x_{1}^{3}x_{2}^{2} - 7680x_{1}^{4}x_{2} + 1344x_{1}^{5} - 4428x_{2}^{6} - 1188x_{1}x_{2}^{5} + 8730x_{1}^{2}x_{2}^{4} + 1240x_{1}^{3}x_{2}^{3} - 5370x_{1}^{4}x_{2}^{2} - 168x_{1}^{5}x_{2} + 952x_{1}^{6} - 648x_{2}^{7} + 1944x_{1}x_{2}^{6} + 3672x_{1}^{2}x_{2}^{5} - 3480x_{1}^{3}x_{2}^{4} - 4080x_{1}^{4}x_{2}^{3} + 2592x_{1}^{5}x_{2}^{2} + 1344x_{1}^{6}x_{2} - 768x_{1}^{7} + 729x_{2}^{8} + 972x_{1}x_{2}^{7} - 1458x_{1}^{2}x_{2}^{6} - 1836x_{1}^{3}x_{2}^{5} + 1305x_{1}^{4}x_{2}^{4} + 1224x_{1}^{5}x_{2}^{3} - 648x_{1}^{6}x_{2}^{2} - 288x_{1}^{7}x_{2} + 144x_{1}^{8} \]
Constraints f(x) - γ to be a sum of squares
con_ref = @constraint(model, f >= γ)
@objective(model, Max, γ)
optimize!(model)-------------------------------------------------------------
Clarabel.jl v0.6.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 45
constraints = 121
nnz(P) = 0
nnz(A) = 121
cones (total) = 2
: Zero = 1, numel = 1
: PSDTriangle = 1, numel = 120
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 6.0000e+02 6.0000e+02 0.00e+00 6.51e-01 4.57e-01 1.00e+00 1.02e+04 ------
1 -3.0355e+03 -2.4055e+03 2.62e-01 1.69e-01 1.33e-01 6.30e+02 2.58e+03 8.06e-01
2 -4.5486e+02 -2.1780e+02 1.09e+00 2.38e-02 2.39e-02 2.37e+02 5.36e+02 8.67e-01
3 3.2891e+01 7.6205e+01 1.32e+00 3.59e-03 3.47e-03 4.33e+01 9.46e+01 8.76e-01
4 2.7894e+01 3.5220e+01 2.63e-01 6.65e-04 6.22e-04 7.33e+00 1.84e+01 8.59e-01
5 1.2781e+01 1.9102e+01 4.95e-01 3.55e-04 4.58e-04 6.32e+00 1.16e+01 5.36e-01
6 5.6656e+00 6.5505e+00 1.56e-01 5.37e-05 7.07e-05 8.85e-01 1.88e+00 8.72e-01
7 4.1810e+00 4.5446e+00 8.70e-02 1.90e-05 2.68e-05 3.64e-01 6.80e-01 8.09e-01
8 3.4403e+00 3.5518e+00 3.24e-02 6.41e-06 8.82e-06 1.12e-01 2.39e-01 7.54e-01
9 3.0921e+00 3.1107e+00 6.02e-03 1.22e-06 1.64e-06 1.86e-02 4.73e-02 9.05e-01
10 3.0531e+00 3.0642e+00 3.64e-03 8.85e-07 9.56e-07 1.11e-02 2.73e-02 5.76e-01
11 3.0065e+00 3.0087e+00 7.22e-04 6.12e-07 1.38e-07 2.17e-03 4.48e-03 9.05e-01
12 3.0037e+00 3.0049e+00 3.99e-04 4.87e-07 7.62e-08 1.20e-03 2.48e-03 6.15e-01
13 3.0009e+00 3.0012e+00 9.46e-05 1.27e-07 1.82e-08 2.84e-04 6.17e-04 8.42e-01
14 3.0003e+00 3.0004e+00 3.35e-05 8.47e-08 6.35e-09 1.01e-04 2.15e-04 7.57e-01
15 3.0000e+00 3.0001e+00 4.74e-06 5.73e-08 8.82e-10 1.42e-05 3.25e-05 9.28e-01
16 3.0000e+00 3.0001e+00 4.74e-06 5.73e-08 8.82e-10 1.42e-05 3.25e-05 0.00e+00
---------------------------------------------------------------------------------------------
Terminated with status = solved (reduced accuracy)
solve time = 28.0msThe lower bound found is 3
solution_summary(model)* Solver : Dual model with Clarabel attached
* Status
Result count : 1
Termination status : ALMOST_OPTIMAL
Message from the solver:
"ALMOST_SOLVED"
* Candidate solution (result #1)
Primal status : NEARLY_FEASIBLE_POINT
Dual status : NEARLY_FEASIBLE_POINT
Objective value : 3.00006e+00
Dual objective value : 3.00004e+00
* Work counters
Solve time (sec) : 2.79524e-02
Barrier iterations : 16
The moment matrix is as follows, we can already see the global minimizer [0, -1] from the entries (2, 1) and (3, 1). This heuristic way to obtain solutions to the polynomial optimization problem is suggested in (Laurent, 2008; (6.15)).
ν = moment_matrix(con_ref)MomentMatrix with row/column basis:
MonomialBasis([1, x[2], x[1], x[2]^2, x[1]*x[2], x[1]^2, x[2]^3, x[1]*x[2]^2, x[1]^2*x[2], x[1]^3, x[2]^4, x[1]*x[2]^3, x[1]^2*x[2]^2, x[1]^3*x[2], x[1]^4])
And entries in a 15×15 SymMatrix{Float64}:
1.0000000018198698 -0.9999904338316166 … 3.4104034285866996e-6
-0.9999904338316166 0.9999813837943897 1.7043839977191587e-6
1.401293336357993e-5 -1.3799592446004206e-5 7.451483754417287e-6
0.9999813860035811 -0.9999722309458006 0.0002711976680818313
-1.3823074641574693e-5 1.4109607254594264e-5 0.00028496583206889375
8.140490999870537e-7 2.3020192377404057e-7 … 0.0005630993629247136
-0.9999722273334571 0.9999633200853798 -0.00013925168205436442
1.4102740426892929e-5 -1.3953042571734061e-5 2.2983476463909564e-6
2.287994538559344e-7 3.914539624501102e-7 0.0005027103259083649
1.14756778339497e-6 9.527302713889533e-7 0.0011791392891342552
0.9999633103973404 -0.999954406161818 … 0.6372092054366799
-1.3906964660954489e-5 1.38684805592237e-5 -0.1551786920868488
3.931301909082682e-7 -2.020302007611492e-7 0.8777398620413609
9.574747243383608e-7 2.4810748752536637e-7 0.2059975266673279
3.4104034285866996e-6 1.7043839977191587e-6 1.422084768180202Many entries of the matrix actually have the same moment. We can obtain the following list of these moments without duplicates (ignoring when difference of entries representing the same moments is below 1e-5)
μ = measure(ν, atol = 1e-5)Measure{Float64, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}([1.0000000018198698, -0.9999904338316166, 1.401293336357993e-5, 0.9999813860035811, -1.3823074641574693e-5, 8.140490999870537e-7, -0.9999722273334571, 1.4102740426892929e-5, 2.287994538559344e-7, 1.14756778339497e-6 … 0.0011791392891342552, 1.4837695259617898, -0.38704799673282014, 0.5307619036335185, -0.3153797748347552, 0.6372092054366799, -0.1551786920868488, 0.8777398620413609, 0.2059975266673279, 1.422084768180202], DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}[1, x₂, x₁, x₂², x₁x₂, x₁², x₂³, x₁x₂², x₁²x₂, x₁³ … x₁⁷, x₂⁸, x₁x₂⁷, x₁²x₂⁶, x₁³x₂⁵, x₁⁴x₂⁴, x₁⁵x₂³, x₁⁶x₂², x₁⁷x₂, x₁⁸])The truncated moment matrix can then be obtained as follows
ν_truncated = moment_matrix(μ, monomials(x, 0:3))MomentMatrix with row/column basis:
MonomialBasis([1, x[2], x[1], x[2]^2, x[1]*x[2], x[1]^2, x[2]^3, x[1]*x[2]^2, x[1]^2*x[2], x[1]^3])
And entries in a 10×10 SymMatrix{Float64}:
1.0000000018198698 -0.9999904338316166 … 1.14756778339497e-6
-0.9999904338316166 0.9999813860035811 9.574747243383608e-7
1.401293336357993e-5 -1.3823074641574693e-5 3.4104034285866996e-6
0.9999813860035811 -0.9999722273334571 2.4810748752536637e-7
-1.3823074641574693e-5 1.4102740426892929e-5 1.7043839977191587e-6
8.140490999870537e-7 2.287994538559344e-7 … 7.451483754417287e-6
-0.9999722273334571 0.9999633103973404 9.773522078850986e-5
1.4102740426892929e-5 -1.3906964660954489e-5 0.0002711976680818313
2.287994538559344e-7 3.931301909082682e-7 0.00028496583206889375
1.14756778339497e-6 9.574747243383608e-7 0.0005630993629247136Let's check if the flatness property is satisfied. The rank of ν_truncated seems to be 1:
using LinearAlgebra
LinearAlgebra.svdvals(Matrix(ν_truncated.Q))
LinearAlgebra.rank(Matrix(ν_truncated.Q), rtol = 1e-3)
svdvals(Matrix(ν_truncated.Q))10-element Vector{Float64}:
3.9999129595584804
0.0008735448885684927
0.00017574267052825453
7.615603475243653e-6
6.052439999404675e-7
3.320433002738225e-7
1.1346221574951914e-7
7.16763637145575e-8
3.140108856202916e-8
4.238028226828164e-9The rank of ν is clearly higher than 1, closer to 3:
svdvals(Matrix(ν.Q))15-element Vector{Float64}:
5.17674761577772
2.3189470566744617
1.453558648851051
0.0023462616989670462
0.0006350998616532194
0.00015172127913365465
0.0001398235480826501
9.002841822263696e-6
6.652718664350959e-7
2.7408971775687256e-7
7.656812574907482e-8
6.435800433409038e-8
1.0665802363009165e-8
6.706830024591651e-9
3.6941439018433997e-10Even if the flatness property is not satisfied, we can still try extracting the minimizer with a low rank decomposition of rank 3. We find the optimal solution again doing so:
atomic_measure(ν, FixedRank(3))Atomic measure on the variables x[1], x[2] with 1 atoms:
at [1.4273511730133842e-5, -0.9999906623763443] with weight 1.032287146482791This page was generated using Literate.jl.